1. Field of the Invention
The present invention relates generally to spacecraft attitude control systems and methods that employ Kalman filters.
2. Description of the Related Art
Kalman filtering is a statistical technique that combines a knowledge of the statistical nature of system measurement errors with a knowledge of system dynamics, as represented in a state space model, to arrive at an estimate of the state of a system. In general, the system state can include any number of unknowns. Kalman filtering for a spacecraft attitude-control system, for example, is typically configured to have a state matrix that includes attitude and gyroscopic bias but much larger dimensions (e.g., on the order of 60) of state matrices are sometimes used (to enhance generality, the term matrices will be used herein even though some matrices may typically be single column or single row matrices that would otherwise be referred to as vectors).
The Kalman filtering process utilizes a weighting function, called the Kalman gain, which is optimized to produce a minimum estimate variance (i.e., the estimate""s accuracy is maximized). In particular, a Kalman filter combines a current measurement y(tn) of a parameter x (e.g., an attitude) at a time tn with measurement and state predictions y*(tnxe2x88x92) and x*(tnxe2x88x92) of the parameter x that are based on past measurements (and thus apply to a time tnxe2x88x92 just before tn) to provide a filtered estimate x*(tn+) of x at a time tn+ just after the time tn. As indicated by the subscript n, the filter successively and recursively combines the measurements and predictions to obtain an estimate with a minimum variance (i.e., maximum accuracy).
This process is succinctly summarized in an estimate update equation
x*(tn+)=x*(tnxe2x88x92)+k(tn){y(tn)xe2x88x92y*(tnxe2x88x92)},xe2x80x83xe2x80x83(1)
in which the state prediction x*(tnxe2x88x92) just before the measurement y(tn) is updated by a portion k(tn) of a residue which is the difference {y(tn)xe2x88x92y*(tnxe2x88x92)} between the measurement y(tn) and the measurement prediction y*(tnxe2x88x92) to form an estimate x*(tn+) for a time tn+ just after the measurement y(tn) was made. The portion k(tn) is the Kalman gain which is calculated as                               k          ⁡                      (                          t              n                        )                          =                                            σ              x                        *                                                           2                            ⁢                              (                                  t                  n                  -                                )                                                                                        σ                x                            *                                                                   2                                ⁢                                  (                                      t                    n                    -                                    )                                                      -                          σ              m              2                                                          (        2        )            
in which "sgr"x*2(tnxe2x88x92) is the estimate variance (i.e., uncertainty of the estimate) just before the measurement y(tn) and "sgr"m2 is the measurement variance (i.e., uncertainty of the measurement). The measurement variance "sgr"m2 is a function of the system under consideration and, more particularly, of the system""s measurement hardware (e.g., manufacturers of spacecraft star trackers and gyros typically specify attitude and attitude rate measurement errors).
In contrast, the estimate variance is reduced as the Kalman process continues. In the beginning of the process, the estimate variance is generally much greater than the measurement variance so that the gain k(tn) of equation (2) approaches one. As the process continues, the estimate variance is reduced below the measurement variance so that the gain k(tn) declines to a value much less than one. It is apparent from equation (1), therefore, that a large portion of the residue {y(tn)xe2x88x92y*(tnxe2x88x92)} is initially used to update the state prediction x*(tnxe2x88x92) into the updated estimate x*(tn+) but this portion decreases as the process continues (i.e., the weighting given to new measurements is successively reduced). In particular, the estimate variance is reduced (i.e., updated) at each measurement y(tn) in accordance with
"sgr"x*2(tn+)=(1xe2x88x92k(tn))"sgr"x*2(tnxe2x88x92)xe2x80x83xe2x80x83(3)
The updated variance "sgr"x*2(tn+) is time delayed so that it becomes the estimate variance "sgr"x*2(tnxe2x88x92) that is used in equation (2) for calculating the next gain k(tn). It has been shown that the estimate variance can be expressed as "sgr"m2/n and thus, it asymptotically approaches zero as more data (i.e., measurements) is obtained.
FIG. 1 illustrates a block diagram 20 of typical Kalman processes that are expressed in a more general matrix form. The Kalman filter 20 comprises an estimator 22 and a gain calculator 24 that supplies a Kalman gain matrix K(tn) to a multiplier 26 of the estimator. To facilitate a description of the filter 20, an investigation of the gain calculator 24 is preceded by the following description of the estimator.
The estimator 22 receives a measurement matrix Y(tn) at an input port 28 and provides an estimate matrix X*(tn+) for a time just after the measurement to an output port 30. From this estimate matrix, a state prediction matrix X*(tnxe2x88x92) is formed for a time just before the next measurement and this state prediction is provided to a summer 32. Because the state of the system typically varies dynamically between measurements, the estimate matrix X*(tn+) that corresponds to a time just after the last measurement must be extrapolated over time to form the state prediction matrix X*(tnxe2x88x92).
As shown in FIG. 1, this extrapolation is accomplished by passing the estimate matrix X*(tn+) through a delay 34 (to cause it to be time-incident with the next measurement matrix Y(tn)) and multiplying it by a state transition matrix 36 which contains extrapolation information in the form of a state transition matrix "PHgr"(tn,tnxe2x88x921). For example, if one component of the estimate matrix X*(tn+) is a position x(tn+), the state transition matrix "PHgr"(tn,tnxe2x88x921) might be configured to extrapolate the position with a term of Tv(tn+) wherein T is the time between measurements and v(tn+) is the last estimate of velocity.
The state prediction matrix X*(tnxe2x88x92) is also multiplied by a measurement matrix 38 to form a measurement prediction Y*(tnxe2x88x92) which is provided to a differencer 40 where it is differenced with the measurement matrix Y(tn). The measurement matrix H(tn) conditions the state prediction matrix X*(tnxe2x88x92) so that its elements correspond to those of the measurement matrix Y(tn) and can be properly differenced with it.
As a first conditioning example, the measurement may be expressed in one coordinate system (e.g., rectangular) and the estimates tracked in a different coordinate system (e.g., spherical). In this example, the measurement matrix H(tn) would be configured to convert the estimates to the coordinate system of the measurements. In a second conditioning example, attitude and gyroscopic bias might be part of the prediction matrix X*(tnxe2x88x92) but only attitude might be present in the measurement matrix Y(tn) so that the measurement matrix H(tn) would be configured to make the necessary conversion.
The differencer 40, therefore, generates a residue Y(tn)xe2x88x92H(tn)Y*(tnxe2x88x92) which is then multiplied in the multiplier 26 to form a correction K(tn){Y(tn)xe2x88x92H(tn)Y(tnxe2x88x92)} that will be used to update the estimate matrix. The updating is performed in the summer 32 where the correction is summed with the state prediction matrix X*(tn) to generate the updated estimate matrix
X*(tn+)=X*(tnxe2x88x92)+K(tn) {Y(tn)xe2x88x92H(tn)Y*(tnxe2x88x92)xe2x80x83xe2x80x83(4)
at the output port 30.
Attention is now directed to the Kalman gain calculator 24 which performs similar updating and extrapolation processes for the Kalman gain matrix K(tn). In the matrix notation of FIG. 1, the estimate variance "sgr"x*2(tnxe2x88x92) and the measurement variance "sgr"m2 are respectively replaced by an estimate covariance matrix P(tn) and a measurement-noise covariance matrix R(tn). These matrices may be respectively considered to be measures of prediction error (i.e., error in the state prediction matrix X*(tnxe2x88x92)) and of measurement error (i.e., error in the estimate matrix Y(tn)).
When these matrix replacements are made in equation (2), the calculation of Kalman gain becomes
K(tn)=P(tnxe2x88x92)HT(tn){H(tn)P(tnxe2x88x92)HT(tn)+R(tn)}xe2x88x92xe2x80x83xe2x80x83(5)
in which the the estimate covariance matrix P(tn) is appropriately modified by the measurement matrix H(tn) and its transform HT(tn). In the gain calculator 24 of FIG. 1, the process recited in equation (5) is performed by the gain calculation block 42.
The estimate covariance matrix P(tn) is updated between measurements in an update block 44 wherein the estimate variance "sgr"x*2(tnxe2x88x92) of equation (3) has been replaced by the estimate covariance matrix P(tn) as modified by the measurement matrix H(tn). Accordingly, the updated estimate covariance matrix P(tn) is processed as
P(tn)={Ixe2x88x92K(tn)H(tn) P(tnxe2x88x92).xe2x80x83xe2x80x83(6)
The updated estimate covariance matrix is passed through delay 45 to cause it to be time-incident with the next measurement matrix Y(tn). Because of the dynamic nature of the system""s state, the updated estimate covariance matrix P(tn) must also be extrapolated over time to form an estimate covariance matrix P(tnxe2x88x92) for the time just before the next measurement. As shown in extrapolation block 46, this extrapolation is processed as
P(tnxe2x88x92)="PHgr"(tn,tnxe2x88x921)P(tn+)"PHgr"T(tn,tnxe2x88x921)+Q(tnxe2x88x921)xe2x80x83xe2x80x83(7)
and is effected with the state transition matrix "PHgr"(tn,tnxe2x88x921) and its transform. In addition, plant or process noise, which includes modeling errors as well as actual noise and system disturbances, is summarized in a process-noise covariance matrix Q(tnxe2x88x921).
As detailed above, the gain calculator 24 updates the estimate covariance matrix P(tn) in update block 44, delays it in delay 45 and extrapolates it in extrapolation block 46 to generate the covariance matrix P(tnxe2x88x92) The updated, delayed and extrapolated estimate covariance matrix in extrapolation block 46 is then combined with the measurement-noise covariance matrix R(tn) in the gain calculation block 42 to calculate a gain K(tn) which is provided to the estimator 22 for processing the current measurement Y(tn) at the input port 28 into an updated estimate X*(tn+) at the output port 30.
Although the Kalman filter 20 of FIG. 1 is configured to combine parameter measurements and parameter predictions to obtain parameter estimates with minimum variances, its operational performance relies upon the accuracy of information that is provided to it (e.g., by the state transition matrix "PHgr"(tn,tnxe2x88x921), the measurement matrix H(tn) and the measurement-noise covariance matrix R(tn)).
In particular, the measurement-noise covariance matrix R(tn) will degrade performance of the Kalman filter of FIG. 1 if it misstates or omits measurement errors because this corrupts the gain K(tn) and, thereby, the calculated estimate X*(tn+). It has been found especially difficult, for example, to accurately model star tracker measurement errors in spacecraft attitude control systems.
The present invention is directed to stellar attitude-control systems and methods whose accuracy is enhanced because they accurately define star tracker measurement variances in Kalman filter gain calculations. These goals are realized with the recognition that important star tracker errors exhibit a boresight symmetry and that these errors can be accurately defined by weighting a measurement-noise covariance matrix R(tn) with variances that are functions of off-boresight angles of detected stars.
A control method embodiment of the invention derives off-boresight angles xcex8 from star-tracker signals of detected stars. These off-boresight angles xcex8 are combined with variance coefficients xcex1 to generate off-boresight variances rob(tn) that are functions of the off-boresight angles xcex8 and, in particular, correspond to star tracker color shift errors.
A gain matrix K(tn) is then calculated with a weighted measurement-noise covariance matrix R(tn) that includes the off-boresight variances rob(tn). This gain matrix is used to generate an attitude estimate matrix X*(tn) for use in spacecraft attitude control.
The variance coefficients xcex1 are preferably modified to include star tracker focal length shift errors which are also a function of off-boresight angles xcex8. It is noted that the variance coefficients a can be reduced by identifying the spectral classes of detected and identified stars.
The invention recognizes that other star tracker errors are not functions of off-boresight angles and therefore modifies the covariance matrix R(tn) with appropriate constant variances rc(tn).
Attitude-control system embodiments are also described for realizing these methods.
The novel features of the invention are set forth with particularity in the appended claims. The invention will be best understood from the following description when read in conjunction with the accompanying drawings.